Data

The data has been collected during the study of . The dataset is available in open-access at …

library(tidyverse)
library(ggside)
library(easystats)
library(patchwork)
library(mirt)

df <- read.csv("https://raw.githubusercontent.com/RealityBending/IllusionGameReliability/main/data/preprocessed_questionnaires.csv") |>
# df <- read.csv("../IllusionGameReliability/data/preprocessed_questionnaires.csv") |> 
  select(Sex, Age, starts_with("Item_PHQ4")) |> 
  filter(!is.na(Item_PHQ4_Anxiety_1))

names(df) <- str_remove_all(names(df), "Item_PHQ4_")

Participants

report::report_participants(df)
## [1] "485 participants (Mean age = 30.1, SD = 10.1, range: [18, 73]; Sex: 50.3% females, 49.7% males, 0.0% other)"

Distribution

add_labels <- function(x) {
  x <- case_when(
      x == 0 ~ "Not at all",
      x == 1 ~ "Once or twice",
      x == 2 ~ "Several days",
      x == 3 ~ "More than half the days",
      TRUE ~ "Nearly every day"
    )
  fct_relevel(x, c("Not at all", "Once or twice", "Several days", "More than half the days",  "Nearly every day"))
}
df <- select(df, -Sex, -Age)
  
data <- df |> 
  pivot_longer(everything(), names_to = "Item", values_to = "Answer") |> 
  group_by(Item, Answer) |> 
  summarise(n = n() / nrow(df)) |> 
  mutate(
    Facet = str_split_fixed(Item, "_", 2)[, 1],
    # Item = str_split_fixed(Item, "_", 2)[,2],
    Answer = add_labels(Answer),
    Item = case_when(
      Item == "Anxiety_1" ~ "A1: Feeling nervous, anxious or on edge",
      Item == "Anxiety_2" ~ "A2: Not being able to stop or control worrying",
      Item == "Depression_3" ~ "D1: Little interest or pleasure in doing things",
      TRUE ~ "D2: Feeling down, depressed, or hopeless"
    ))

p1 <- data |> 
  ggplot(aes(x = Answer, fill=Facet)) +
  geom_bar(aes(y = n), stat = "identity") +
  scale_y_continuous(labels=scales::percent, expand = c(0, 0)) +
  scale_fill_manual(values = c("Depression" = "#3F51B5", "Anxiety" = "#E91E63"), guide = "none") +
  labs(y = "Proportion of Answers", title = "A. Proportion of answers per item") +
  theme_modern(axis.title.space = 10) +
  theme(axis.text.x = element_text(angle = 45, hjust=1),
        axis.title.x = element_blank(),
        plot.title = element_text(face = "bold"),
        strip.background=element_rect(fill="#EEEEEE", colour="white")) +
  facet_wrap(~Item)
p1

The “Once or twice” answer was selected in 29.85%

p1b <- rbind(
  data.frame(Scores = paste0(df$Anxiety_1, "_", df$Anxiety_2),
             Facet = "Anxiety"),
  data.frame(Scores = paste0(df$Depression_3, "_", df$Depression_4),
             Facet = "Depression")
) |> 
  group_by(Facet, Scores) |> 
  summarize(n = n() / nrow(df)) |> 
  separate(Scores, into = c("Q1", "Q2")) |> 
  mutate(Label = ifelse(Q1 == Q2, format_percent(n), "")) |> 
  mutate(Q1 = add_labels(as.numeric(Q1)),
         Q2 = add_labels(as.numeric(Q2))) |> 
  ggplot(aes(x = Q1, y = Q2)) +
  geom_tile(aes(fill = n)) +
  geom_text(aes(label = Label)) +
  scale_y_discrete(expand = c(0, 0)) +
  scale_x_discrete(expand = c(0, 0)) +
  scale_fill_gradientn(colors=c("white", "#2196F3", "#3F51B5", "#673AB7"), labels = scales::percent) +
  labs(fill = "% Participants", x = "Item 1", y = "Item 2", title = "Joint prevalence of responses") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust=1),
        panel.grid.major = element_blank(),
        plot.title = element_text(face = "bold"),
        strip.background=element_rect(fill="#EEEEEE", colour="white"),
        strip.text = element_text(face = "bold")) +
  facet_wrap(~Facet) 
p1b

Anxiety

Descriptive

df_anx <- select(df, contains("Anxiety"))

performance::cronbachs_alpha(df_anx)
## [1] 0.903

Model

m_anxiety <- mirt(df_anx, model = 1, itemtype = 'graded', SE = TRUE)
## 
Iteration: 1, Log-Lik: -1434.510, Max-Change: 0.49722
Iteration: 2, Log-Lik: -1377.689, Max-Change: 0.62182
Iteration: 3, Log-Lik: -1322.367, Max-Change: 0.72396
Iteration: 4, Log-Lik: -1283.817, Max-Change: 0.66639
Iteration: 5, Log-Lik: -1263.174, Max-Change: 0.50226
Iteration: 6, Log-Lik: -1253.813, Max-Change: 0.34052
Iteration: 7, Log-Lik: -1249.772, Max-Change: 0.23422
Iteration: 8, Log-Lik: -1247.939, Max-Change: 0.16157
Iteration: 9, Log-Lik: -1247.027, Max-Change: 0.11472
Iteration: 10, Log-Lik: -1246.280, Max-Change: 0.06322
Iteration: 11, Log-Lik: -1246.106, Max-Change: 0.05105
Iteration: 12, Log-Lik: -1246.000, Max-Change: 0.04207
Iteration: 13, Log-Lik: -1245.837, Max-Change: 0.01752
Iteration: 14, Log-Lik: -1245.829, Max-Change: 0.01703
Iteration: 15, Log-Lik: -1245.823, Max-Change: 0.01663
Iteration: 16, Log-Lik: -1245.792, Max-Change: 0.01528
Iteration: 17, Log-Lik: -1245.788, Max-Change: 0.01454
Iteration: 18, Log-Lik: -1245.783, Max-Change: 0.01500
Iteration: 19, Log-Lik: -1245.760, Max-Change: 0.01522
Iteration: 20, Log-Lik: -1245.755, Max-Change: 0.01616
Iteration: 21, Log-Lik: -1245.750, Max-Change: 0.01580
Iteration: 22, Log-Lik: -1245.726, Max-Change: 0.01663
Iteration: 23, Log-Lik: -1245.722, Max-Change: 0.01803
Iteration: 24, Log-Lik: -1245.718, Max-Change: 0.01645
Iteration: 25, Log-Lik: -1245.694, Max-Change: 0.02301
Iteration: 26, Log-Lik: -1245.689, Max-Change: 0.01862
Iteration: 27, Log-Lik: -1245.684, Max-Change: 0.01712
Iteration: 28, Log-Lik: -1245.660, Max-Change: 0.02176
Iteration: 29, Log-Lik: -1245.656, Max-Change: 0.02120
Iteration: 30, Log-Lik: -1245.651, Max-Change: 0.02021
Iteration: 31, Log-Lik: -1245.626, Max-Change: 0.01921
Iteration: 32, Log-Lik: -1245.622, Max-Change: 0.02052
Iteration: 33, Log-Lik: -1245.618, Max-Change: 0.01939
Iteration: 34, Log-Lik: -1245.594, Max-Change: 0.02217
Iteration: 35, Log-Lik: -1245.589, Max-Change: 0.02285
Iteration: 36, Log-Lik: -1245.585, Max-Change: 0.02256
Iteration: 37, Log-Lik: -1245.561, Max-Change: 0.01765
Iteration: 38, Log-Lik: -1245.557, Max-Change: 0.01948
Iteration: 39, Log-Lik: -1245.553, Max-Change: 0.02142
Iteration: 40, Log-Lik: -1245.529, Max-Change: 0.02326
Iteration: 41, Log-Lik: -1245.525, Max-Change: 0.02321
Iteration: 42, Log-Lik: -1245.521, Max-Change: 0.02237
Iteration: 43, Log-Lik: -1245.499, Max-Change: 0.02324
Iteration: 44, Log-Lik: -1245.495, Max-Change: 0.02257
Iteration: 45, Log-Lik: -1245.491, Max-Change: 0.02463
Iteration: 46, Log-Lik: -1245.468, Max-Change: 0.02322
Iteration: 47, Log-Lik: -1245.464, Max-Change: 0.02381
Iteration: 48, Log-Lik: -1245.461, Max-Change: 0.02434
Iteration: 49, Log-Lik: -1245.440, Max-Change: 0.02352
Iteration: 50, Log-Lik: -1245.436, Max-Change: 0.02455
Iteration: 51, Log-Lik: -1245.432, Max-Change: 0.02511
Iteration: 52, Log-Lik: -1245.411, Max-Change: 0.02529
Iteration: 53, Log-Lik: -1245.408, Max-Change: 0.02560
Iteration: 54, Log-Lik: -1245.405, Max-Change: 0.02584
Iteration: 55, Log-Lik: -1245.385, Max-Change: 0.02519
Iteration: 56, Log-Lik: -1245.381, Max-Change: 0.02795
Iteration: 57, Log-Lik: -1245.378, Max-Change: 0.02268
Iteration: 58, Log-Lik: -1245.365, Max-Change: 0.02836
Iteration: 59, Log-Lik: -1245.361, Max-Change: 0.02557
Iteration: 60, Log-Lik: -1245.358, Max-Change: 0.02714
Iteration: 61, Log-Lik: -1245.340, Max-Change: 0.02686
Iteration: 62, Log-Lik: -1245.337, Max-Change: 0.02687
Iteration: 63, Log-Lik: -1245.334, Max-Change: 0.02683
Iteration: 64, Log-Lik: -1245.316, Max-Change: 0.02909
Iteration: 65, Log-Lik: -1245.313, Max-Change: 0.02611
Iteration: 66, Log-Lik: -1245.311, Max-Change: 0.02433
Iteration: 67, Log-Lik: -1245.296, Max-Change: 0.02794
Iteration: 68, Log-Lik: -1245.293, Max-Change: 0.02868
Iteration: 69, Log-Lik: -1245.290, Max-Change: 0.02824
Iteration: 70, Log-Lik: -1245.274, Max-Change: 0.03027
Iteration: 71, Log-Lik: -1245.272, Max-Change: 0.03022
Iteration: 72, Log-Lik: -1245.269, Max-Change: 0.02880
Iteration: 73, Log-Lik: -1245.253, Max-Change: 0.01696
Iteration: 74, Log-Lik: -1245.252, Max-Change: 0.03171
Iteration: 75, Log-Lik: -1245.249, Max-Change: 0.02984
Iteration: 76, Log-Lik: -1245.238, Max-Change: 0.02126
Iteration: 77, Log-Lik: -1245.234, Max-Change: 0.02745
Iteration: 78, Log-Lik: -1245.231, Max-Change: 0.00321
Iteration: 79, Log-Lik: -1245.231, Max-Change: 0.02698
Iteration: 80, Log-Lik: -1245.229, Max-Change: 0.02876
Iteration: 81, Log-Lik: -1245.226, Max-Change: 0.02913
Iteration: 82, Log-Lik: -1245.213, Max-Change: 0.02962
Iteration: 83, Log-Lik: -1245.211, Max-Change: 0.02906
Iteration: 84, Log-Lik: -1245.209, Max-Change: 0.02921
Iteration: 85, Log-Lik: -1245.196, Max-Change: 0.03105
Iteration: 86, Log-Lik: -1245.194, Max-Change: 0.02980
Iteration: 87, Log-Lik: -1245.192, Max-Change: 0.03064
Iteration: 88, Log-Lik: -1245.179, Max-Change: 0.02993
Iteration: 89, Log-Lik: -1245.177, Max-Change: 0.03042
Iteration: 90, Log-Lik: -1245.175, Max-Change: 0.03243
Iteration: 91, Log-Lik: -1245.164, Max-Change: 0.02865
Iteration: 92, Log-Lik: -1245.161, Max-Change: 0.03030
Iteration: 93, Log-Lik: -1245.159, Max-Change: 0.03026
Iteration: 94, Log-Lik: -1245.149, Max-Change: 0.02627
Iteration: 95, Log-Lik: -1245.147, Max-Change: 0.00253
Iteration: 96, Log-Lik: -1245.147, Max-Change: 0.03315
Iteration: 97, Log-Lik: -1245.145, Max-Change: 0.02951
Iteration: 98, Log-Lik: -1245.144, Max-Change: 0.03432
Iteration: 99, Log-Lik: -1245.142, Max-Change: 0.03095
Iteration: 100, Log-Lik: -1245.134, Max-Change: 0.02948
Iteration: 101, Log-Lik: -1245.132, Max-Change: 0.03010
Iteration: 102, Log-Lik: -1245.131, Max-Change: 0.03003
Iteration: 103, Log-Lik: -1245.122, Max-Change: 0.02659
Iteration: 104, Log-Lik: -1245.120, Max-Change: 0.02996
Iteration: 105, Log-Lik: -1245.118, Max-Change: 0.03053
Iteration: 106, Log-Lik: -1245.110, Max-Change: 0.03140
Iteration: 107, Log-Lik: -1245.108, Max-Change: 0.03240
Iteration: 108, Log-Lik: -1245.107, Max-Change: 0.03022
Iteration: 109, Log-Lik: -1245.099, Max-Change: 0.03116
Iteration: 110, Log-Lik: -1245.097, Max-Change: 0.02909
Iteration: 111, Log-Lik: -1245.096, Max-Change: 0.03552
Iteration: 112, Log-Lik: -1245.091, Max-Change: 0.02542
Iteration: 113, Log-Lik: -1245.087, Max-Change: 0.02774
Iteration: 114, Log-Lik: -1245.085, Max-Change: 0.00788
Iteration: 115, Log-Lik: -1245.085, Max-Change: 0.03110
Iteration: 116, Log-Lik: -1245.084, Max-Change: 0.03182
Iteration: 117, Log-Lik: -1245.082, Max-Change: 0.03100
Iteration: 118, Log-Lik: -1245.075, Max-Change: 0.02973
Iteration: 119, Log-Lik: -1245.074, Max-Change: 0.00171
Iteration: 120, Log-Lik: -1245.074, Max-Change: 0.03305
Iteration: 121, Log-Lik: -1245.073, Max-Change: 0.03095
Iteration: 122, Log-Lik: -1245.072, Max-Change: 0.03123
Iteration: 123, Log-Lik: -1245.070, Max-Change: 0.02556
Iteration: 124, Log-Lik: -1245.066, Max-Change: 0.03212
Iteration: 125, Log-Lik: -1245.064, Max-Change: 0.03086
Iteration: 126, Log-Lik: -1245.063, Max-Change: 0.02970
Iteration: 127, Log-Lik: -1245.057, Max-Change: 0.02988
Iteration: 128, Log-Lik: -1245.056, Max-Change: 0.02452
Iteration: 129, Log-Lik: -1245.055, Max-Change: 0.02659
Iteration: 130, Log-Lik: -1245.050, Max-Change: 0.03050
Iteration: 131, Log-Lik: -1245.049, Max-Change: 0.00127
Iteration: 132, Log-Lik: -1245.049, Max-Change: 0.02366
Iteration: 133, Log-Lik: -1245.048, Max-Change: 0.03117
Iteration: 134, Log-Lik: -1245.047, Max-Change: 0.03019
Iteration: 135, Log-Lik: -1245.046, Max-Change: 0.02695
Iteration: 136, Log-Lik: -1245.042, Max-Change: 0.02767
Iteration: 137, Log-Lik: -1245.041, Max-Change: 0.02447
Iteration: 138, Log-Lik: -1245.040, Max-Change: 0.02816
Iteration: 139, Log-Lik: -1245.036, Max-Change: 0.02278
Iteration: 140, Log-Lik: -1245.035, Max-Change: 0.00219
Iteration: 141, Log-Lik: -1245.035, Max-Change: 0.03010
Iteration: 142, Log-Lik: -1245.034, Max-Change: 0.02924
Iteration: 143, Log-Lik: -1245.033, Max-Change: 0.02585
Iteration: 144, Log-Lik: -1245.032, Max-Change: 0.02333
Iteration: 145, Log-Lik: -1245.029, Max-Change: 0.03156
Iteration: 146, Log-Lik: -1245.028, Max-Change: 0.02978
Iteration: 147, Log-Lik: -1245.028, Max-Change: 0.00251
Iteration: 148, Log-Lik: -1245.028, Max-Change: 0.00250
Iteration: 149, Log-Lik: -1245.028, Max-Change: 0.03614
Iteration: 150, Log-Lik: -1245.027, Max-Change: 0.03004
Iteration: 151, Log-Lik: -1245.023, Max-Change: 0.00262
Iteration: 152, Log-Lik: -1245.023, Max-Change: 0.02886
Iteration: 153, Log-Lik: -1245.022, Max-Change: 0.02489
Iteration: 154, Log-Lik: -1245.019, Max-Change: 0.03097
Iteration: 155, Log-Lik: -1245.018, Max-Change: 0.00207
Iteration: 156, Log-Lik: -1245.018, Max-Change: 0.00074
Iteration: 157, Log-Lik: -1245.018, Max-Change: 0.00062
Iteration: 158, Log-Lik: -1245.018, Max-Change: 0.00066
Iteration: 159, Log-Lik: -1245.018, Max-Change: 0.00015
Iteration: 160, Log-Lik: -1245.018, Max-Change: 0.00013
Iteration: 161, Log-Lik: -1245.018, Max-Change: 0.00049
Iteration: 162, Log-Lik: -1245.018, Max-Change: 0.00059
Iteration: 163, Log-Lik: -1245.018, Max-Change: 0.00019
Iteration: 164, Log-Lik: -1245.018, Max-Change: 0.00040
Iteration: 165, Log-Lik: -1245.018, Max-Change: 0.00014
Iteration: 166, Log-Lik: -1245.018, Max-Change: 0.00012
Iteration: 167, Log-Lik: -1245.018, Max-Change: 0.00048
Iteration: 168, Log-Lik: -1245.018, Max-Change: 0.00010
## 
## Calculating information matrix...

m_anxiety
## 
## Call:
## mirt(data = df_anx, model = 1, itemtype = "graded", SE = TRUE)
## 
## Full-information item factor analysis with 1 factor(s).
## Converged within 1e-04 tolerance after 168 EM iterations.
## mirt version: 1.37.1 
## M-step optimizer: BFGS 
## EM acceleration: Ramsay 
## Number of rectangular quadrature: 61
## Latent density type: Gaussian 
## 
## Information matrix estimated with method: Oakes
## Second-order test: model is a possible local maximum
## Condition number of information matrix =  30371
## 
## Log-likelihood = -1245
## Estimated parameters: 10 
## AIC = 2510
## BIC = 2552; SABIC = 2520
## G2 (14) = 20.6, p = 0.113
## RMSEA = 0.031, CFI = NaN, TLI = NaN

# Factor loadings
summary(m_anxiety)
##              F1    h2
## Anxiety_1 0.895 0.802
## Anxiety_2 0.991 0.982
## 
## SS loadings:  1.78 
## Proportion Var:  0.892 
## 
## Factor correlations: 
## 
##    F1
## F1  1

# Alpha
a <- coef(m_anxiety)
a <- data.frame(
  Item = c("Anxiety 1", "Anxiety 2"),
  Discrimination = c(paste0(format_value(a$Anxiety_1[1, 1]), ", ", format_ci(a$Anxiety_1[2, 1], a$Anxiety_1[3, 1])), paste0(format_value(a$Anxiety_2[1, 1]), ", ", format_ci(a$Anxiety_2[2, 1], a$Anxiety_2[3, 1]))
))
knitr::kable(a)
Item Discrimination
Anxiety 1 3.42, 95% CI [2.66, 4.18]
Anxiety 2 12.55, 95% CI [-9.55, 34.65]
# Plots
plot(m_anxiety, type = 'trace', theta_lim = c(-3, 3))
plot(m_anxiety, type = 'infotrace', theta_lim = c(-3, 3))
plot(m_anxiety, type = 'score', theta_lim = c(-3, 3))
plot(m_anxiety, type = 'infoSE', theta_lim = c(-3, 3))
plot(m_anxiety, type = 'rxx', theta_lim = c(-3, 3))

Category Characteristic Curves (CRC) and Item Information Curves

Typically, an optimally informative item will have a large location and broad category coverage (as indicated by location parameters) over theta.

crc <- function(mod, data) {
  Theta <- matrix(seq(-3.5, 3.5, by = .01))
  rez <- data.frame()
  for(i in 1:2) {
    dat <- as.data.frame(probtrace(extract.item(mod, i), Theta))
    dat$Theta <- Theta[,1]
    dat$Information <- normalize(iteminfo(extract.item(mod, i), Theta))
    dat <- pivot_longer(dat, -one_of(c("Theta", "Information")), names_to = "Answer", values_to = "Probability")
    dat$Item <- names(data)[i]
    rez <- rbind(rez, dat)
  }

  rez <- rez |>
    mutate(
      Answer = case_when(
      Answer == "P.1" ~ "Not at all",
      Answer == "P.2" ~ "Once or twice",
      Answer == "P.3" ~ "Several days",
      Answer == "P.4" ~ "More than half the days",
      TRUE ~ "Nearly every day"
      ),
      Answer = fct_relevel(Answer, c("Not at all", "Once or twice", "Several days", "More than half the days",  "Nearly every day")),
    Item = case_when(
      Item == "Anxiety_1" ~ "A1: Feeling nervous, anxious or on edge",
      Item == "Anxiety_2" ~ "A2: Not being able to stop or control worrying",
      Item == "Depression_3" ~ "D1: Little interest or pleasure in doing things",
      TRUE ~ "D2: Feeling down, depressed, or hopeless"
    ))
  
  # Get close to 0.95
  i1 <- rez[rez$Item == unique(rez$Item)[1],]
  i2 <- rez[rez$Item == unique(rez$Item)[2],]
  minmax <- rbind(
    rbind(
      i1[i1$Answer == "Not at all", ][which.min(abs(i1[i1$Answer == "Not at all", ]$Probability - 0.95)), ],
      i1[i1$Answer == "Nearly every day", ][which.min(abs(i1[i1$Answer == "Nearly every day", ]$Probability - 0.95)), ]
      ),
    rbind(
      i2[i2$Answer == "Not at all", ][which.min(abs(i2[i2$Answer == "Not at all", ]$Probability - 0.95)), ],
      i2[i2$Answer == "Nearly every day", ][which.min(abs(i2[i2$Answer == "Nearly every day", ]$Probability - 0.95)), ]
      )
  )

  p <- rez |>
    ggplot(aes(x = Theta, y = Probability)) +
    geom_line(aes(y = Information), linetype = "dotted", color = "grey") +
    geom_line(aes(color = Answer), size=1) +
    # geom_vline(data = minmax, aes(xintercept = Theta), linetype = "dashed") +
    scale_y_continuous(labels=scales::percent, expand = c(0, 0.005)) +
    scale_color_flat_d("rainbow") +
    facet_grid(~Item) +
    theme_modern(axis.title.space = 5) +
    theme(strip.background=element_rect(fill="#EEEEEE", colour="white"))
  list(p = p, rez = rez, minmax=minmax)
}

out <- crc(m_anxiety, df_anx)
p2a <- out$p + labs(x = expression(Anxiety~(theta)))
p2a

Normalized scoring

normalize_scores <- function(data, minmax) {

  minmax <- minmax[minmax$Theta %in% range(minmax$Theta), ]
  
  item <- data.frame()
  scores <- data.frame()
  for(i in unique(data$Item)) {
    item <- data[data$Item == i, ] |>
      filter(Theta >= min(minmax$Theta) & Theta <= max(minmax$Theta)) |>
      mutate(Theta_n = normalize(Theta)) |> 
      rbind(item)

    scores <- item[item$Item == i, ] |>
      group_by(Answer) |>
      filter(Probability == max(Probability)) |>
      ungroup() |>
      mutate(label = insight::format_value(Theta_n)) |> 
      rbind(scores)
  }
  
  item |>
    ggplot(aes(x = Theta, y = Probability, color = Answer)) +
    geom_line(size = 1) +
    geom_segment(data=scores, aes(x=Theta, xend = Theta, y=0, yend=Probability, color = Answer), linetype = "dashed") +
    geom_label(data=scores, aes(x=Theta, y=Probability, label=label), show.legend = FALSE) +
    scale_y_continuous(labels=scales::percent, expand = c(0, 0.01)) +
    scale_color_viridis_d(option = "D") +
    facet_grid(~Item) +
    theme_modern(axis.title.space = 5) +
    theme(strip.background=element_rect(fill="#EEEEEE", colour="white"))
}

p3a <- normalize_scores(data=out$rez, out$minmax) + labs(x = expression(Anxiety~(theta)))
p3a 

Depression

Descriptive

df_dep <- select(df, contains("Depression"))

performance::cronbachs_alpha(df_dep)
## [1] 0.841

Model

m_depression <- mirt(df_dep, model = 1, itemtype = 'graded', SE = TRUE)
## 
Iteration: 1, Log-Lik: -1454.429, Max-Change: 0.43021
Iteration: 2, Log-Lik: -1412.776, Max-Change: 0.48366
Iteration: 3, Log-Lik: -1375.753, Max-Change: 0.50038
Iteration: 4, Log-Lik: -1351.861, Max-Change: 0.41059
Iteration: 5, Log-Lik: -1339.730, Max-Change: 0.29809
Iteration: 6, Log-Lik: -1334.341, Max-Change: 0.20880
Iteration: 7, Log-Lik: -1332.019, Max-Change: 0.14321
Iteration: 8, Log-Lik: -1330.986, Max-Change: 0.09907
Iteration: 9, Log-Lik: -1330.503, Max-Change: 0.07020
Iteration: 10, Log-Lik: -1330.136, Max-Change: 0.03563
Iteration: 11, Log-Lik: -1330.079, Max-Change: 0.02845
Iteration: 12, Log-Lik: -1330.048, Max-Change: 0.02329
Iteration: 13, Log-Lik: -1330.002, Max-Change: 0.01230
Iteration: 14, Log-Lik: -1329.997, Max-Change: 0.01214
Iteration: 15, Log-Lik: -1329.991, Max-Change: 0.01163
Iteration: 16, Log-Lik: -1329.957, Max-Change: 0.01344
Iteration: 17, Log-Lik: -1329.950, Max-Change: 0.01509
Iteration: 18, Log-Lik: -1329.943, Max-Change: 0.01385
Iteration: 19, Log-Lik: -1329.900, Max-Change: 0.02264
Iteration: 20, Log-Lik: -1329.890, Max-Change: 0.01837
Iteration: 21, Log-Lik: -1329.882, Max-Change: 0.01814
Iteration: 22, Log-Lik: -1329.825, Max-Change: 0.02072
Iteration: 23, Log-Lik: -1329.815, Max-Change: 0.02112
Iteration: 24, Log-Lik: -1329.804, Max-Change: 0.02158
Iteration: 25, Log-Lik: -1329.734, Max-Change: 0.02587
Iteration: 26, Log-Lik: -1329.720, Max-Change: 0.02565
Iteration: 27, Log-Lik: -1329.707, Max-Change: 0.02592
Iteration: 28, Log-Lik: -1329.623, Max-Change: 0.03023
Iteration: 29, Log-Lik: -1329.607, Max-Change: 0.03013
Iteration: 30, Log-Lik: -1329.591, Max-Change: 0.03048
Iteration: 31, Log-Lik: -1329.494, Max-Change: 0.03517
Iteration: 32, Log-Lik: -1329.477, Max-Change: 0.03493
Iteration: 33, Log-Lik: -1329.459, Max-Change: 0.01506
Iteration: 34, Log-Lik: -1329.456, Max-Change: 0.04136
Iteration: 35, Log-Lik: -1329.437, Max-Change: 0.03616
Iteration: 36, Log-Lik: -1329.419, Max-Change: 0.03661
Iteration: 37, Log-Lik: -1329.309, Max-Change: 0.04138
Iteration: 38, Log-Lik: -1329.290, Max-Change: 0.04095
Iteration: 39, Log-Lik: -1329.271, Max-Change: 0.04120
Iteration: 40, Log-Lik: -1329.157, Max-Change: 0.04555
Iteration: 41, Log-Lik: -1329.137, Max-Change: 0.04525
Iteration: 42, Log-Lik: -1329.117, Max-Change: 0.04549
Iteration: 43, Log-Lik: -1329.002, Max-Change: 0.04943
Iteration: 44, Log-Lik: -1328.983, Max-Change: 0.04984
Iteration: 45, Log-Lik: -1328.964, Max-Change: 0.04835
Iteration: 46, Log-Lik: -1328.853, Max-Change: 0.05603
Iteration: 47, Log-Lik: -1328.834, Max-Change: 0.05338
Iteration: 48, Log-Lik: -1328.816, Max-Change: 0.05192
Iteration: 49, Log-Lik: -1328.712, Max-Change: 0.05430
Iteration: 50, Log-Lik: -1328.694, Max-Change: 0.05488
Iteration: 51, Log-Lik: -1328.677, Max-Change: 0.05551
Iteration: 52, Log-Lik: -1328.578, Max-Change: 0.05994
Iteration: 53, Log-Lik: -1328.562, Max-Change: 0.05838
Iteration: 54, Log-Lik: -1328.546, Max-Change: 0.05845
Iteration: 55, Log-Lik: -1328.456, Max-Change: 0.06069
Iteration: 56, Log-Lik: -1328.441, Max-Change: 0.06055
Iteration: 57, Log-Lik: -1328.427, Max-Change: 0.06088
Iteration: 58, Log-Lik: -1328.345, Max-Change: 0.06152
Iteration: 59, Log-Lik: -1328.332, Max-Change: 0.06212
Iteration: 60, Log-Lik: -1328.319, Max-Change: 0.06246
Iteration: 61, Log-Lik: -1328.245, Max-Change: 0.06481
Iteration: 62, Log-Lik: -1328.233, Max-Change: 0.06452
Iteration: 63, Log-Lik: -1328.222, Max-Change: 0.06456
Iteration: 64, Log-Lik: -1328.155, Max-Change: 0.06593
Iteration: 65, Log-Lik: -1328.145, Max-Change: 0.06591
Iteration: 66, Log-Lik: -1328.134, Max-Change: 0.06874
Iteration: 67, Log-Lik: -1328.073, Max-Change: 0.06981
Iteration: 68, Log-Lik: -1328.063, Max-Change: 0.06728
Iteration: 69, Log-Lik: -1328.054, Max-Change: 0.06611
Iteration: 70, Log-Lik: -1328.003, Max-Change: 0.07361
Iteration: 71, Log-Lik: -1327.994, Max-Change: 0.06790
Iteration: 72, Log-Lik: -1327.986, Max-Change: 0.06682
Iteration: 73, Log-Lik: -1327.940, Max-Change: 0.06446
Iteration: 74, Log-Lik: -1327.933, Max-Change: 0.07094
Iteration: 75, Log-Lik: -1327.926, Max-Change: 0.07007
Iteration: 76, Log-Lik: -1327.884, Max-Change: 0.07049
Iteration: 77, Log-Lik: -1327.877, Max-Change: 0.07142
Iteration: 78, Log-Lik: -1327.871, Max-Change: 0.08124
Iteration: 79, Log-Lik: -1327.831, Max-Change: 0.06560
Iteration: 80, Log-Lik: -1327.824, Max-Change: 0.06915
Iteration: 81, Log-Lik: -1327.818, Max-Change: 0.07046
Iteration: 82, Log-Lik: -1327.786, Max-Change: 0.07482
Iteration: 83, Log-Lik: -1327.780, Max-Change: 0.07337
Iteration: 84, Log-Lik: -1327.775, Max-Change: 0.07155
Iteration: 85, Log-Lik: -1327.746, Max-Change: 0.07769
Iteration: 86, Log-Lik: -1327.741, Max-Change: 0.07462
Iteration: 87, Log-Lik: -1327.736, Max-Change: 0.07441
Iteration: 88, Log-Lik: -1327.710, Max-Change: 0.07091
Iteration: 89, Log-Lik: -1327.706, Max-Change: 0.07165
Iteration: 90, Log-Lik: -1327.702, Max-Change: 0.07258
Iteration: 91, Log-Lik: -1327.678, Max-Change: 0.07570
Iteration: 92, Log-Lik: -1327.675, Max-Change: 0.07507
Iteration: 93, Log-Lik: -1327.671, Max-Change: 0.07409
Iteration: 94, Log-Lik: -1327.651, Max-Change: 0.07329
Iteration: 95, Log-Lik: -1327.647, Max-Change: 0.07495
Iteration: 96, Log-Lik: -1327.644, Max-Change: 0.07514
Iteration: 97, Log-Lik: -1327.625, Max-Change: 0.07639
Iteration: 98, Log-Lik: -1327.622, Max-Change: 0.07650
Iteration: 99, Log-Lik: -1327.619, Max-Change: 0.00312
Iteration: 100, Log-Lik: -1327.619, Max-Change: 0.07501
Iteration: 101, Log-Lik: -1327.616, Max-Change: 0.07486
Iteration: 102, Log-Lik: -1327.613, Max-Change: 0.07507
Iteration: 103, Log-Lik: -1327.597, Max-Change: 0.07353
Iteration: 104, Log-Lik: -1327.594, Max-Change: 0.07149
Iteration: 105, Log-Lik: -1327.592, Max-Change: 0.00414
Iteration: 106, Log-Lik: -1327.592, Max-Change: 0.00415
Iteration: 107, Log-Lik: -1327.591, Max-Change: 0.00021
Iteration: 108, Log-Lik: -1327.591, Max-Change: 0.07427
Iteration: 109, Log-Lik: -1327.589, Max-Change: 0.07333
Iteration: 110, Log-Lik: -1327.587, Max-Change: 0.07407
Iteration: 111, Log-Lik: -1327.584, Max-Change: 0.07417
Iteration: 112, Log-Lik: -1327.570, Max-Change: 0.07582
Iteration: 113, Log-Lik: -1327.568, Max-Change: 0.07539
Iteration: 114, Log-Lik: -1327.566, Max-Change: 0.05673
Iteration: 115, Log-Lik: -1327.561, Max-Change: 0.00505
Iteration: 116, Log-Lik: -1327.560, Max-Change: 0.08820
Iteration: 117, Log-Lik: -1327.558, Max-Change: 0.07629
Iteration: 118, Log-Lik: -1327.546, Max-Change: 0.00593
Iteration: 119, Log-Lik: -1327.545, Max-Change: 0.00122
Iteration: 120, Log-Lik: -1327.545, Max-Change: 0.07510
Iteration: 121, Log-Lik: -1327.544, Max-Change: 0.00080
Iteration: 122, Log-Lik: -1327.544, Max-Change: 0.07048
Iteration: 123, Log-Lik: -1327.542, Max-Change: 0.06928
Iteration: 124, Log-Lik: -1327.532, Max-Change: 0.06730
Iteration: 125, Log-Lik: -1327.530, Max-Change: 0.00117
Iteration: 126, Log-Lik: -1327.530, Max-Change: 0.00300
Iteration: 127, Log-Lik: -1327.530, Max-Change: 0.00188
Iteration: 128, Log-Lik: -1327.530, Max-Change: 0.00069
Iteration: 129, Log-Lik: -1327.530, Max-Change: 0.00076
Iteration: 130, Log-Lik: -1327.530, Max-Change: 0.00018
Iteration: 131, Log-Lik: -1327.530, Max-Change: 0.00055
Iteration: 132, Log-Lik: -1327.530, Max-Change: 0.00057
Iteration: 133, Log-Lik: -1327.530, Max-Change: 0.00023
Iteration: 134, Log-Lik: -1327.530, Max-Change: 0.00048
Iteration: 135, Log-Lik: -1327.530, Max-Change: 0.00014
Iteration: 136, Log-Lik: -1327.530, Max-Change: 0.00011
Iteration: 137, Log-Lik: -1327.530, Max-Change: 0.00251
Iteration: 138, Log-Lik: -1327.530, Max-Change: 0.00075
Iteration: 139, Log-Lik: -1327.530, Max-Change: 0.00094
Iteration: 140, Log-Lik: -1327.530, Max-Change: 0.00015
Iteration: 141, Log-Lik: -1327.530, Max-Change: 0.00045
Iteration: 142, Log-Lik: -1327.530, Max-Change: 0.00058
Iteration: 143, Log-Lik: -1327.530, Max-Change: 0.00016
Iteration: 144, Log-Lik: -1327.530, Max-Change: 0.00043
Iteration: 145, Log-Lik: -1327.530, Max-Change: 0.00093
Iteration: 146, Log-Lik: -1327.530, Max-Change: 0.00054
Iteration: 147, Log-Lik: -1327.530, Max-Change: 0.00027
Iteration: 148, Log-Lik: -1327.530, Max-Change: 0.00022
Iteration: 149, Log-Lik: -1327.530, Max-Change: 0.00038
Iteration: 150, Log-Lik: -1327.530, Max-Change: 0.00013
Iteration: 151, Log-Lik: -1327.530, Max-Change: 0.00011
Iteration: 152, Log-Lik: -1327.530, Max-Change: 0.00041
Iteration: 153, Log-Lik: -1327.530, Max-Change: 0.00038
Iteration: 154, Log-Lik: -1327.530, Max-Change: 0.00021
Iteration: 155, Log-Lik: -1327.530, Max-Change: 0.00041
Iteration: 156, Log-Lik: -1327.530, Max-Change: 0.00013
Iteration: 157, Log-Lik: -1327.530, Max-Change: 0.00011
Iteration: 158, Log-Lik: -1327.530, Max-Change: 0.00188
Iteration: 159, Log-Lik: -1327.530, Max-Change: 0.00071
Iteration: 160, Log-Lik: -1327.530, Max-Change: 0.00278
Iteration: 161, Log-Lik: -1327.530, Max-Change: 0.00053
Iteration: 162, Log-Lik: -1327.530, Max-Change: 0.00011
Iteration: 163, Log-Lik: -1327.530, Max-Change: 0.00009
## 
## Calculating information matrix...

m_depression
## 
## Call:
## mirt(data = df_dep, model = 1, itemtype = "graded", SE = TRUE)
## 
## Full-information item factor analysis with 1 factor(s).
## Converged within 1e-04 tolerance after 163 EM iterations.
## mirt version: 1.37.1 
## M-step optimizer: BFGS 
## EM acceleration: Ramsay 
## Number of rectangular quadrature: 61
## Latent density type: Gaussian 
## 
## Information matrix estimated with method: Oakes
## Second-order test: model is a possible local maximum
## Condition number of information matrix =  72063
## 
## Log-likelihood = -1328
## Estimated parameters: 10 
## AIC = 2675
## BIC = 2717; SABIC = 2685
## G2 (14) = 27.1, p = 0.0188
## RMSEA = 0.044, CFI = NaN, TLI = NaN

# Factor loadings
summary(m_depression)
##                 F1    h2
## Depression_3 0.995 0.989
## Depression_4 0.817 0.667
## 
## SS loadings:  1.66 
## Proportion Var:  0.828 
## 
## Factor correlations: 
## 
##    F1
## F1  1

# Alpha
a2 <- coef(m_depression)
a2 <- data.frame(
  Item = c("Depression 1", "Depression 2"),
  Discrimination = c(paste0(format_value(a2$Depression_3[1, 1]), ", ", format_ci(a2$Depression_3[2, 1], a2$Depression_3[3, 1])), paste0(format_value(a2$Depression_4[1, 1]), ", ", format_ci(a2$Depression_4[2, 1], a2$Depression_4[3, 1]))
))
knitr::kable(a2)
Item Discrimination
Depression 1 16.46, 95% CI [-10.08, 43.00]
Depression 2 2.41, 95% CI [2.03, 2.79]

Category Characteristic Curves (CRC) and Item Information Curves

out <- crc(m_depression, df_dep)
p2b <- out$p + labs(x = expression(Depression~(theta)))
p2b

Normalized scoring

p3b <- normalize_scores(out$rez, out$minmax) + labs(x = expression(Depression~(theta)))
p3b 

Figures

fig1 <- wrap_elements(p1) + wrap_elements(p2a / p2b + plot_annotation(title = "B. Information Curves per Response Type", theme = list(plot.title = element_text(face = "bold"))) + plot_layout(guides = "collect")) + plot_layout(widths = c(1, 1.5))

ggsave("figures/figure1.png", fig1, width = fig.width*1.43, height=fig.height*1.43)

fig2 <- p3a / p3b + plot_annotation(title = "Normalized Scoring", theme = list(plot.title = element_text(face = "bold"))) + plot_layout(guides = "collect")

ggsave("figures/figure2.png", fig2, width = fig.width*1.43, height=fig.height*1.43)

Scores

data <- sapply(cbind(df_anx, df_dep), function(x) {
  case_when(
      x == 0 ~ "Not at all",
      x == 1 ~ "Once or twice",
      x == 2 ~ "Several days",
      x == 3 ~ "More than half the days",
      TRUE ~ "Nearly every day"
    )
}) |> 
  as.data.frame()

source("score_PHQ4.R")

data <- score_PHQ4(A1=data$Anxiety_1, A2=data$Anxiety_2, D1=data$Depression_3, D2=data$Depression_4, method="basic") |> 
  datawizard::data_addsuffix("_Basic") |> 
  cbind(
    score_PHQ4(A1=data$Anxiety_1, A2=data$Anxiety_2, D1=data$Depression_3, D2=data$Depression_4, method="normalized") |> 
      datawizard::data_addsuffix("_Normalized")) |> 
  select(matches("Depression|Anxiety")) |> 
  pivot_longer(everything()) |> 
  separate(name, into = c("Dimension", "Scoring")) 

data |> 
  ggplot(aes(x=value)) + 
  geom_area(data=estimate_density(data, at = c("Dimension", "Scoring"), method = "KernSmooth") |> 
              group_by(Scoring) |> 
              normalize(select = "y"),
            aes(x = x, y = y, fill=Dimension), alpha = 0.05, position="identity") +
  geom_hline(yintercept = 0.5, linetype = "dotted") +
  stat_ecdf(aes(color=Dimension), geom = "step", size=1) +
  facet_wrap(~Scoring, scales = "free") +
  scale_y_continuous(labels=scales::percent, expand = c(0, 0)) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_fill_manual(values = c("Depression" = "#3F51B5", "Anxiety" = "#E91E63"), guide = "none") +
  scale_color_manual(values = c("Depression" = "#3F51B5", "Anxiety" = "#E91E63")) +
  theme_modern(axis.title.space = 5) +
  labs(title = "Distribution and Cumulative Density of the Sample", x = "Score", y = "Cumulative Proportion")

Correlation

df$Depression <- mirt::fscores(m_depression)[, 1]
df$Anxiety <- mirt::fscores(m_anxiety)[, 1]

correlation::correlation(df) |>
  summary() |> 
  export_table(caption = "Correlation Matrix")
## Correlation Matrix
## 
## Parameter    | Anxiety | Depression | Depression_4 | Depression_3 | Anxiety_2
## -----------------------------------------------------------------------------
## Anxiety_1    |    0.89 |       0.72 |         0.56 |         0.72 |      0.83
## Anxiety_2    |    0.98 |       0.72 |         0.57 |         0.73 |          
## Depression_3 |    0.74 |       0.98 |         0.73 |              |          
## Depression_4 |    0.59 |       0.78 |              |              |          
## Depression   |    0.75 |            |              |              |